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Abstract. Synchrony is one of the most common dynamical states emerging on networks. The speed of con- 
vergence towards synchrony provides a fundamental collective time scale for synchronizing systems. Here 
we study the asymptotic synchronization times for directed networks with topologies ranging from com- 
pletely ordered, grid-like, to completely disordered, random, including intermediate, partially disordered 
topologies. We extend the approach of Master Stability Functions to quantify synchronization times. We 
find that the synchronization times strongly and systematically depend on the network topology. In partic- 
ular, at fixed in-degree, stronger topological randomness induces faster synchronization, whereas at fixed 
path length, synchronization is slowest for intermediate randomness in the small-world regime. Randomly 
rewiring real-world neural, social and transport networks confirms this picture. 
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1 Introduction and overview 

We live in a world where everything is connected: we are 
surrounded by global networks of communication, trans- 
portation, trade, social relations and media. In addition, 
due to the enormous progress in technology, we are now 
able to decipher natural networks of high complexity. For 
example, we get to know more details about the brain 
or the human genome day by day. But understanding 
the function of these networks requires a two-sided ap- 
proach. Firstly, we need to know the basic structures of 
these networks at both microscopic and macroscopic level. 
Secondly, we have to map the rules which govern the dy- 
namic interactions. Although all these networks seem to 
be different at first sight, there are similarities such as ab- 
stract patterns or simple organizing principles [T] . It is not 
only the individual unit that matters but also the archi- 
tecture of the connections. To uncover these rules, theo- 
retical studies first focus on systems consisting of simple 
units such as oscillators and of simple structures: all-to-all 
coupled units, lattices, or mean field models. Although ex- 
act results on synchronization in networks with a general 
structure have been obtained recently [2 , 3 , 4 J , it is still not 
well understood how the structure of a complex network 
affects dynamical features of synchronization. 

But a question which dates back to the sixties changed 
this approach, namely: 'what is the probability that any 
two people, selected arbitrarily from a large population, 
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such as that of the United States, will know each other?' 
|S]. A more interesting formulation, however, takes into 
account that, while two persons may not know each other 
directly, they may share one or more mutual acquaintances 
such that any two people are connected through a chain 
of acquaintances. This concept, known as 'six degrees of 
separation', refers to the idea that everyone is on aver- 
age approximately six steps away from any other person 
on Earth. In 1998, this idea was converted into a simple 
model of small- world networks [BJ. The crucial point is 
that this model interpolates between totally regular and 
totally random topologies. 

Starting with a ring where units only communicate 
with their direct neighbours, ring connections between 
neighbours are cut and connected to randomly chosen 
nodes somewhere else in the ring. An important quantity 
which characterizes this network's evolving architecture is 
the 'average path length' and formalizes the intuitive idea 
of degrees of separation. Due to the rewiring, the path 
length between any pair of units drastically decreases. 
Simultaneously this architecture still exhibits high local 
clustering, defined as the probability that two nodes linked 
to a common node will also be linked to each other. These 
two properties were suggested to be particularly support- 
ive of synchronization. Indeed, several detailed studies sup- 
port this view by showing that at fixed coupling strength 
small-world networks tend to already synchronize at lower 
connectivity than many other classes of networks 0E] . 
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Synchronization is one of the most frequently observed 
collective dynamics in many physical and biological sys- 
tems [TJ[71[S]. Synchronization might be both advantageous 
and desired, for instance in secure communication |9J, or 
detrimental and undesired, as during tremor in patients 
with Parkinson disease or during epileptic seizures |101 

El. 

Therefore, a broad area of research has emerged [12; 13 
I14| , determining under which conditions on the interaction 
strengths and topologies coupled units actually synchro- 
nize and when they do not. These results suggest some 
key properties about the topological influence on the net- 
work synchronizability, i.e. the capability of a network to 
synchronize at all. They do not tell much about the speed 
of synchronization given that a network synchronizes in 
principle. 

For any real system, however, it equally matters how 
fast the units synchronize or whether the network interac- 
tions fail to coordinate the units' dynamics on time scales 
relevant to the system's function (or dysfunction) , cf . |15l 
IT51fT7lfT5| . The applications range from consensus dynam- 
ics of distributed decision-making problems for interacting 
groups of agents [TU] to neuroscience questions of how fast 
the visual processing or olfactory discrimination could be 
[20 , 21] . Yet this question is far from being understood and 
currently under active investigation [22 , 23 . 24 , 25 . 26ll2"7] . 
In particular it is largely unknown how fast small worlds 
synchronize which leads us to our main question addressed 
here: what is the typical time scale for synchronization, 
i.e. how fast can oscillators coordinate their dynamics if 
they are not directly interconnected but interact on large 
networks of regular, random or small- world topology? 

We address this question by computer simulations as 
well as analytical predictions. All results are derived for 
the simplest of all regular states, the synchronous periodic 
state, in which all oscillators exhibit identical dynamics. 
However, also other settings are imaginable: cluster states 
in which two or more groups of synchronized oscillators 
exist [28 , 29 or systems with inhomogeneities in the dy- 
namical and topological parameters |30| can be treated 
similarly. 

We study the effect of topology on the synchronization 
time of directed networks which exhibit different dynam- 
ics: Kuramoto phase oscillators coupled via phase differ- 
ences, higher-dimensional periodic Rossler systems cou- 
pled diffusively as well as neural circuits with inhibitory 
delayed pulse-coupling. Synchronization time is a mea- 
sure of how quickly the network synchronizes after being 
perturbed from a synchronized state. So far it has been 
studied analytically for fully random networks only [24] . 
Firstly, comparing network ensembles with a fixed number 
of edges, it is shown that those in the small- world regime 
synchronize faster than regular networks but slower than 
random networks. This is expected intuitively - the char- 
acteristic path length is monotonically decreasing while 
rewiring - and in accordance with the result for synchro- 
nizability Hence, we fix the average characteristic 
path length and again investigate the dependence of syn- 
chronization time on the network's topology. We find that 



- for a fixed average characteristic path length - networks 
in the small-world regime again synchronize slower than 
random networks, but this time even slower than regu- 
lar networks: we see a non-monotonic dependence on the 
topological randomness. First results have been reported 
in [22 and here we further systematically investigate these 
studies and extend them to real-world networks. We com- 
pare network ensembles with fixed topological quantities 
like the betweenness centrality as well as generic ensem- 
bles for Kuramoto, Rossler and pulse-coupled oscillators. 
Moreover, we make analytical predictions of the synchro- 
nization times for periodic Rossler systems and observe 
remarkable similarities between the synchronization times 
for the Kuramoto and pulse-coupled oscillators. 

This article is organized as follows. In Section|2]we first 
introduce the concept of synchronization time, the central 
quantity of our studies. In Section [3] we explain the un- 
derlying network structure and introduce the small-world 
model adapted to directed networks. Section [4] outlines 
the different types of considered dynamics. Here we de- 
rive analytical predictions for the synchronization times 
and extend the master stability function formalism [1"4"] 
to determine the synchronization speed. In Section [5] we 
compare the analytical predictions for the synchronization 
times and the results obtained by extensive computer sim- 
ulations for network ensembles with fixed in-degree, with 
fixed average path length and with fixed betweenness cen- 
trality, followed by an analysis of generic network ensem- 
bles. In Section [6] the study of synchronization times for 
real-world networks - rewiring them towards fully ran- 
dom networks - confirms our theoretical results. We close 
in Section [7] with a summary and a discussion of further 
work. 



2 Synchronization time 

We assume identical oscillators which means we get com- 
plete synchronization in the end. The equation of motion 
for the uncoupled oscillators looks like 

f -*•<*<), a) 

where the m-dimensional vector Xi = {xi i, Xi >m } refers 
to the components of each oscillator i £ {1,...,N} and 
F : R m i-> R m defines the dynamics. 

We describe the connection of N oscillators in a net- 
work by a coupling matrix J that basically consists of zero 
and non-zero elements to specify which oscillators are cou- 
pled to which other ones. The matrix elements are thus 
given by 

j _ j a /ki if j is connected to i ^ j 
13 ' ~ \ otherwise ^ ' 

where a is a global coupling constant and 1/fcj a normal- 
ization factor that guarantees a homogeneous total input 

N 
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Fig. 1. (Color online) Time scales of synchronization of oscillator networks for topological randomness q £ {0, 0.02, 1} (in-degree 
fixed at k — 20). top row a: Kuramoto oscillators (a — 1); b: diffusively coupled periodic Rossler oscillators (a = 0.2, c = 5.7, 
b = 1.7, a — 2); c: pulse-coupled oscillators (/ = 1.01, 7 = 1, a = —0.2, A = 0.1). See equations Q - (111 for the definitions 
of the variable differences. Plotted in the right-hand column are the logarithmized decaying distances (see equations Q, (181, 
(22) and pT)) 



such that no specific oscillator receives distinguished cou- 
plings. 

Directly related to the coupling matrix J is the graph 
Laplacian A denned as 

Aij = J lj {l - <%) - adij , (4) 

where Sij is the Kronecker-delta. Considering directed net- 
works its eigenvalues A^ are complex and ordered as = 
Re Ai > Re A2 > ... > Re \n- The number of zero eigen- 
values of the Laplacian matrix is equal to the number of 
strongly-connected components (SCCs) of the networlQ It 
is trivial then to conclude that if Re A2 = 0, the network 
is split in more than one SCC. Then, from a dynamical 
point of view, it is impossible for the network to achieve 
a complete synchronized state, which is only possible for 
subnetworks with internal coherence. We are considering 
networks with one SCC only throughout this article, which 
means we have always Re A2 > 0. Note that for the pulse- 
coupled system in Section |4~3] the eigenvalues are ordered 
according to their absolute values. 

We describe the dynamics of the interaction with a 
function H that is a vector function of dimension m of 
the dynamical variables of two connected oscillators. Each 
oscillator has the same interaction function. For example, 
H for the Rossler oscillators is a 3 x 3 matrix that only 

1 A strongly-connected component is the maximal subnet- 
work such that for each given pair of nodes (i,j) there is a 
directed path from i to j and from j to i. 



picks out the rr-component to couple to the other oscilla- 
tors. The coupled equations of motion become 

^ = F( Xi ) + JijH(xi, xj) , (5) 

where > acts on each oscillator as a whole. Further- 
more, note again that we only consider identical oscilla- 
tors which means that the global coupling a is the same 
for each oscillator as well as the interaction function H . 
Since we want to examine the case of identical synchro- 
nization, the equations of motion become the same for all 
oscillators when the system is synchronized. In the syn- 
chronous state all oscillators' variables are equal to the 
same dynamical variable: 

X!(t) = x 2 (t) = . . . = x N (t) = s(t) , (6) 

where s(t) is a solution of (|l) as long as H(s(t), s(t)) = 0, 
which is the case for Kuramoto and Rossler oscillators. 
The subspace defined by the constraint of setting all os- 
cillator vectors to the same, synchronous vector is called 
the synchronization manifold. We assume stability of this 
state which means that small arbitrary perturbations to 
each Xj die out in the long time limit . 

In addition to these dynamical systems with continuous- 
time coupling we introduce pulse-coupled systems as well 
in Section |4~3] 

We consider directed regular, small- world and random 
networks which are characterized by increasing rewiring, 
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the topological randomness q. By tuning this parameter 
we interpolate between regular ring networks (q — 0), 
small worlds (low q <C 1) and fully random networks 
(q = 1) which is explained in detail in Section [3] 

First simulations for three different kinds ofoscillators 
(see Fig. [T]) show that synchronization becomes an expo- 
nential process after some short transients for all fractions 
q £ [0, 1] of randomness. Thus the distance 



q = 



q =0.1 



q = 1 



d(t) = m&xdist(xi(t),Xj(t)) 



hJ 



from the synchronous state decays as 
d{t) ~ exp(-t/r) 



(7) 



(8) 



in the long time limit, where dist(a;, x') is a function mea- 
suring the distance between the two appropriate phase 
variables x and x' ', taking into account the periodic bound- 
ary conditions. The characteristic time scale r in pH) is 
what we call the synchronization time in the following. 
Note that there exist systems as well where the transient 
until the exponential decay is not negligible |15 l fT7 j [3T]. 

As one can see in Fig. [I] this decay is similar for Ku- 
ramoto, Rossler and pulse-coupled oscillators. It depicts 
the differences of the phase variables (which we defined in 
detail in Section [4} of ten randomly chosen oscillators to 
the corresponding means denoted by [ . ] : 



A K , i {t) = e i {t)-[e j {t)] 3 



lPC,i 



(t) = &(t) 



,-(*) 



with 



i(t) 



4>i{t) if ^(t) < 0.5 
<pi(t)-l if 4>i{t) > 0.5 



(9) 
(10) 

(11) 
(12) 



'K' stands for Kuramoto, 'R' for Rossler and 'PC for 
pulse-coupled oscillators and these abbreviations will be 
kept throughout this paper. 

In contrast to the continuous-time dynamics of the Ku- 
ramoto and Rossler oscillators, for the pulse-coupled oscil- 
lators the phases are measured at discrete 'spiking' times 
of a reference oscillator. For the 3-dimensional Rossler os- 
cillators only the x-coordinates are shown here. The actual 
dynamical variables for all systems will be introduced in 
Section g] 



3 Network topology 

We adapt the standard small-world model of Watts and 
Strogatz [B] to directed networks |32| . We start with regu- 
lar ring networks where each unit i receives directed links 
from its ki/2 nearest neighbors on both sides (fc is chosen 
to be even). Here the in-degree ki = k Vi is the same for 
all units. We randomly cut each outgoing edge with prob- 
ability q and rewire it to a node chosen uniformly at ran- 
dom from the whole network (avoiding double edges and 
self-loops). We do, however, allow the edge to be rewired 




Fig. 2. (Color online) Rewiring directed networks, the in- 
degree k for each node stays fixed: this can be observed here 
for two reference nodes' incoming edges depicted in blue and 
red. (Here: TV = 10, k = 4, q € {0, 0.1, 1}). 



back to its original position. An important observation 
here is that as q varies the in-degree of each node (and 
with it the average in-degree of the network) is still k (see 
Fig. [2]). This is due to the fact we only rewire outgoing 
egdes. Furthermore the entries > of the coupling 
matrix are multiplied by a global coupling constant a and 
normalized to guarantee that each oscillator i is getting 
the same input as has already been mentioned in Section 
2] The matrix elements are therefore Jjy = er/fci if 



there is a connection from j to i 7^ j, Jij = if there is 
no connection and Ju — for the diagonal elements 

To analyze the purely topological impact on the syn- 
chronization times, we study the network dynamics in its 
simplest setting: we consider strongly-connected networks 
with fixed in-degree k and homogeneous total input cou- 
pling strengths (encoded in the coupling matrix J 
such that full synchrony is achieved from sufficiently close 
initial conditions for all coupling strengths a [25 . 

The directed small-world networks behave as in the 
original Watts-Strogatz model. The small-world regime (a 
small) is characterized by a large clustering coefficienlr] 
(C(q, k)) and a small average path lengtlj^] (L(q, fc)). Here 
( . ) denotes averaging over network realizations at given 
q and fc. 

To quantitatively fix the small world regime we take 



(L(q,k)) 
L(0,k) 



< 0.5 and 



(C(q,k)) 
C(0,fc) 



> 0.85 



(13) 



throughout this study. The results below are not sensitive 
to a change of these values. 

As the topological randomness q is changed from to 
1 the network interpolates between regular and random 
topologies. This structural change induces changes in the 
corresponding graph Laplacian's spectrum and thus has 
a direct influence on the synchronization speed as is ex- 
plained in detail in Section [4j 

2 C(q, k) denotes the actual divided by the possible number 
of directed triangles containing a given node i, averaged over 
all i. 

3 L(q, k) denotes the length of the shortest directed path 
between a given pair of nodes averaged over all 
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4 Oscillator dynamics on networks 

We consider various oscillator types, intrinsic dynamics 
and coupling schemes: phase oscillators coupled via phase 
differences, neural circuits with inhibitory delayed pulse- 
coupling and higher-dimensional periodic systems coupled 
diffusively. In the three following subsections we introduce 
the theory of these different types of oscillators and add 
remarks on the simulations, the chosen initial conditions 
and the numerical measurement of the synchronization 
time. 



A2 here determines the asymptotic synchronization time 
which is given by 

This feature was recently shown to hold also more gener- 
ally for network systems where the stability matrix is not 
necessarily proportional to the graph Laplacian [1,23,34 . 

4.2 Rossler oscillators 



4.1 Kuramoto oscillators 

Consider N Kuramoto oscillators [33 that interact on a 
directed network. Here the dynamical variable of each os- 
cillator is Xi := &i € S 1 = 27rR/N, i.e. a one-dimensional 
phase, with its interaction function H(9i,9j) :— sin(#j — 
6i). Therefore the dynamics of phases 0j(t) of oscillators i 
with time t satisfy 



dt 



%) forze {!,... ,N} , (14) 



where uj is the natural frequency of the oscillators. The 
fully synchronous state defined in ^ here takes the form 



i(t) = 9j(t) =:6{t) 



(15) 



As the synchronous periodic orbit analyzed is isolated 
in state space, the relaxation time continuously changes 
with possible inhomogeneities, so the qualitative results 
obtained below are generic and also hold in the presence 
of small heterogeneities, cf. [30 . 

Furthermore, starting from random initial phases in 
the range [0, 7r] the synchronization dynamics shows a fast 
transient. After this fast initial evolution all phases are 
quite similar and the sine function in ( 14 1 can be well 
approximated by its argument. Linearizing (14 1 close to 



the synchronous state (151 phase perturbations defined as 



$K,i(f) :=0<(t)-0(*) 
evolve according to 



d& 



^=X>i* KlJ -(t) far* € {!,...,#}. 



dt 



(16) 



(17) 



Here the stability matrix coincides with the graph Lapla- 
cian defined in Q. 

Close to every invariant trajectory the eigenvalue A2 
of the stability matrix A that is second largest in real part 
dominates the asymptotic decay in the long time limit 



4(f) ~ exp(-t/T K ) 



(18) 



The distance d^(t) is given by ^ where dist(0, 9') for 
Kuramoto oscillators is the circular distance between the 
two phases 9 and 9' on S 1 . 



Consider a network of Rossler oscillators, both in the chaotic 
and in the periodic regime. Each elementary oscillator is 
described now by three variables {x(i), y(t), z(t)}. The col- 
lective dynamics of N coupled, identical Rossler oscillators 
(i € {1, 2, N}) is governed by the equations 



N 



y, 



3 = 1 



yi =Xi + ayi, 

Zi =b + Zi(xi - c), 



(20) 



where a, b and c are fixed parameters. 

To study the Rossler system in the periodic regime we 
set the parameters to a — 0.2, b = 1.7, c = 5.7. Analo- 
gously setting the parameters toa = 0.2, 6 = 0. 2, c = 5.7 
the chaotic attractor is gained. 

The evolution of perturbations is characterized by mea- 
suring the Euclidean distances 

dy(t) = ^{xm-xm 2 +^{t)-vm 2 +(^{t)-zm 2 

(21) 

between the states of all N(N — l)/2 possible pairs of 
oscillators The asymptotic synchronization time is 

then determined via the decay of the maximal distance 



^r(^) = Taaxdij(t). 



(22) 



A general approach to determine the synchronization 
time for continuous systems described by (JHJ- alternative 
to the one taken for the Kuramoto oscillators, which does 
not work for the Rossler oscillators - is to extend the 
master stability function (MSF) formalism introduced in 
|T4] . Note that this approach does not work for the pulse- 
coupled oscillators, where the phases are measured at dis- 
crete times. So far this formalism has only been used to 
determine the stability of networks of coupled oscillators 
[35 ,36_ and nearly all studies have focussed on symmetric 
undirected networks (see [37J for an exception). 

Defining infinitesimal perturbations to the synchronous 
state (|6j) in the system described by Eq. ^ as 



<*R,i = Xi(t) - S(t) 



(23) 



we get the variational equation 

N 



dS R ,i 
dt 



DF(s)8 R!i -J2A ij DH(s,s)8 Ki , (24) 
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-0.5 



real part a 




Fig. 3. (Color online) Master stability functions (30 1 for the 
Rossler oscillators, a: periodic regime (parameters set to a — 
0.2, b — 1.7, c = 5.7). b: chaotic regime (parameters set to a — 
0.2, b — 0.2, c = 5.7; the black contour shows the MSF equal 
to zero, i.e. separates the stable from the unstable region). 



where the matrix A is the graph Laplacian defined in 
Q, DF(s) and DH(s,s) are the Jacobians evaluated 
along the trajectory s(t). 

For the above Rossler system with diffusive coupling 
via the ^-coordinate the Jacobian matrices for this block 
are given by 



DF(x,y,z) 




and 



DH{x 



1 0> 

//■ i) = | 
,0 0, 



(25) 



(26) 



The transformation S'r = _1 <$r, where O is a matrix 
whose columns are the set of the Laplacian's eigenvectors, 
diagonalizes the set of equations ( 24 1 and hence leads to 



a set of decoupled blocks of the form 



dS' 



R.i 



dt 



[DF(s)-XiDH(s,s)]S' 



(27) 



with the Xi being the eigenvalues of the Laplacian matrix 
A. The above-given Jacobians evaluated in the synchro- 
nized state s(t) are the same for each block, hence the 
blocks only differ by the scalar multiplier A^. 



53. 

t: 
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Fig. 4. (Color online) The graph Laplacian's eigenvalue dis- 
tribution directly links to the system's stability and synchro- 
nization dynamics, a: in gray the whole complex spectrum of 
the graph Laplacians in the ranges q £ [0, 1] and k G [10, 100], 
blue: eigenvalues (purely real due to the initial ring symmetry) 
for q — 0, k = 50, red: eigenvalues in the SW regime with 
q — 0.02, k — 50, green: eigenvalues for q = 1, k = 50. b: zoom 
to the real axis, eigenvalues for k — 50 and q G {0, 0.02, 1} with 
colors as above. 



Thus these blocks could be evaluated all at once by 
setting 



dS'i 



dt 



[DF{s) -{a + 6(3)DH{s, «)] S' R , 



(28) 



in dependence on the complex coupling parameter a + Q/3. 
The imaginary part j3 may be interpreted as a 'rotation' 
taking place between the several decaying eigenmodes of 
the system |38| . 

The system actually synchronizes if 



hi,i = lim - log — - — < 

t^oo t \0 R,i(0)| 



(29) 



for all i € {2, . . . , N}. Here is the largest Lyapunov 
exponent corresponding to the mode of eigenvalue Xi (see 
e.g. E2|). 

To obtain the asymptotic synchronization time we ex- 
tract the largest Lyapunov exponent h\ i with the mini- 
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mal absolute value out of the N 
namely 

hi. 



1 maximal exponents. 

(30) 



= max hi j. 

i>2 

This is also called the master stability function and plotted 
in Fig. [3] for the periodic and chaotic Rossler oscillators. 
We calculate the largest Lyapunov exponent following the 
numerical procedure described in |36| . It is evident that 
the Ai = 0-mode is parallel to the synchronization mani- 
fold while all the other modes are transverse to it. The syn- 
chronization time for the Rossler oscillators is then given 
by 

tr = -j-^- (31) 

"4, max 

^l.max dominates the decay towards the synchronized 
state, but note the nonlinear dependence on the eigen- 
values of the Laplacian matrix (27). Only for the simple 



1-dimcnsional Kuramoto oscillators there is a direct rela- 
tion, since here there is linear and unbounded coupling, 
i.e. the larger the global coupling the faster the synchro- 
nization speed. 

In order to find a value for the global coupling pa- 
rameter a - encoded in A according to Q - that leads 
to synchronization, one calculates the whole spectrum of 
possible eigenvalues (see Fig. [i} to guarantee that each 
one is located in the stable region, i.e. in the region where 
the MSF takes only negative values. Note that the MSF 
for uncoupled periodic Rossler oscillators is zero, while 
the MSF for uncoupled chaotic Rossler oscillators is pos- 
itive. This means that the minimal global coupling con- 
stant needed to achieve synchronization is always larger 
for the chaotic Rossler oscillators than for the periodic 
ones. 



4.3 Pulse-coupled oscillators 

Moreover, we investigated the collective dynamics of pulse- 
coupled neural oscillators [ 171140] , which do not exactly fit 
the general description in Section [2] In this case the dy- 
namical oscillator variables are the membrane potentials 
Vi(t) and delayed discrete output pulses satisfying 



~dt 



N 



Here, each potential Vj relaxes towards I > 1 and is reset 
to zero whenever it reaches a threshold at unity, 



1 



Vj(t) :=0, tj 



t, and to i — ^ to + 1. 



(33) 

At these times tj tm neuron j sends a pulse that after a de- 
lay A > changes the potential of post-synaptic neurons i 
in an inhibitory (negative) manner according to ( 32 1 with 
a < in ((2). 

Equivalent to these ordinary differential equations there 
is a simplified approach which represents the state of a 
one-dimensional oscillator not by its membrane potential, 
but by a phase that encodes the time to the next spike in 



the absence of any interactions. The state of an individ- 
ual oscillator j is then represented by a phase-like variable 
4>j G ( — oo, 1] that increases uniformly in time, 



bj/dt = 1 



(34) 



Upon crossing the firing threshold, 4>j(t~) = 1, at time t 
an oscillator is instantaneously reset to zero, cj)j{t) = 0, 
and a pulse is sent. After a delay time A this pulse is 
received by all oscillators i connected to j and induces an 
instantaneous phase jump given by 

<f>i((t + A)+) = l/- 1 (U(4>i(t + A) + Jij) 

Here, the coupling strengths from j to i are taken to be 
purely inhibitory (er < in ([2])) and normalized according 
to (J3). The rise function U, which mediates the interac- 
tions, can be derived from (32 1 (41], and turns out to be 
monotonic increasing, U' > 0, concave (down), U" < 0, 
and represents the subthreshold dynamics of individual 
oscillators. Note that the function U need to be defined 
on the entire range of accessible phase values. In particu- 
lar, inhibitory coupling can lead to negative phase values 
&<0. ' 

The synchronous state s(t) defined in ^ here takes 
the form 

(35) 

for all i, which is a self-consistent solution assuming that 
all neuronal oscillators fire at the same time. Here all os- 
cillators display identical phases <j>o(t) on a periodic orbit 
such that 4>q (t + T) = <j>o(i) with the period 



where 



T = A+ 1 - a 



= U- 1 (U(A)+a). 



(36) 



(37) 



Note that here in contrast to the Kuramoto oscillators the 
period is different from the one of a free oscillator |3] . 
A perturbation 



£pc(0) =: <5pc = (Sp C ,i, ■ ■ ■ 
to the phases is defined as 

5p C ,i = &(o) - <Mo) 



,Spo,n) (38) 
(39) 



The initial condition for the phases of the pulse-coupled 
oscillators is a random perturbation (5pc from the globally 
synchronized state £pc = 0. The perturbation's compo- 
nents <5pc,i are each drawn independently from a uniform 
distribution on [—6, 5]. The condition 8 < y derived in (24] 
(recall that A is the delay time) ensures that the globally 
synchronized state is stable. This guarantees that all the 
neurons fire before any spikes are received. 

A sufficiently small perturbation 8pc asymptotically 
converges exponentially with time to a constant vector. 
Subtracting the asymptotic phase shift, 



^pcW : = <Spc(*) " Hm S PC (s), 



(40) 
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the distance 



and 



dp C (riT) := max \6' PC i(nT)\ 



from the synchronous state (S' PC i = 0) decays as 

TlT 

dpc(nT) ~ exp( ) 

tpc 



(41) 



(42) 



as n — > oo, defining a synchronization time rpc- 

To understand how the speed of synchronization de- 
pends on the dynamical and network parameters, we an- 
alyze how perturbations 8pc to the synchronous state 
evolve in time. Following [4 we first define a nonlinear 
stroboscopic map 



6 PC (nT) = G(6pc((n-l)T)) 



(43) 



for the perturbations. Note that dpc(T) — £pc(0) since 
no spikes are received before all the oscillators reach the 
phase threshold for the first time. Hence we first apply the 
map G in the first period when spikes are received i.e. for 
n > 2. 

Considering the first order approximation of this period- 
T map one gets a linear iterative map A given by 

N 

SpcAnT) = AjS PC ,j((n - 1)T) , n > 2, (44) 
i=i 

for the perturbations 5pc,i(nT) of spike times close to the 
synchronous orbit of period T = In (//(/ — 1)). 
The matrix elements An are defined as 



Pi m—i if j is connected to i ^ j 

if j = i " (45) 

otherwise 



where the variables p, iTO (to S {1, . . . , ki}) encode phase 
jumps evoked by all pulses up to the mth one received 
jl]. Since the matrix elements (45| are differences of these 
Pi t m , matrix elements Aij and ~A\ji with j ^ j' have in 
general different values depending on the order of incom- 
ing signals. 

This multi-operator problem [32] is induced by the 
structure of the network together with the pulsed inter- 
actions, in particular, by the order of the components of 
6(0). For networks with homogeneous, global coupling dif- 
ferent matrices A can be identified by an appropriate per- 
mutation of the oscillator indices. But in general this is 
impossible. However, here we focus on the integrate-and- 
fire dynamics where the matrix A becomes independent 
of the rank order of the perturbations [22] . Here U takes 
the form 




U(cf>) := -(1 
7 



-i<t> 



(46) 



In order to obtain the matrix elements A^ we first calcu- 
late 



U'(<j>) = le'^ 



(47) 



U-\y) = tl n (l-yi)-i 

7 l 



7 



(48) 



Furthermore we calculate 



U- 1 (U(A) + J ij ) = i]n(e-^ - p^)" 1 (49) 

and 

This leads to 



U'{U-\U{A) + Jij)) = Ie-~< A - 7J tf . (50) 



Pi 



U>(U-i(U(A) + Jij)) 
1 c 7 2^,1=1 J jji 



(51) 



where the sum ^i. J i j l with me {1, . . . , ki} counts up to 
the mth signal received by neuron i during the considered 
period. 

For homogeneous inhibitory coupling, a/k < for each 
existing connection, the elements of the stability matrix 
are given by 

do jk if j is connected to i ^ j 
A i:j = { l-a iij = i (52) 
otherwise 



with 



a 



Ie~~< A + 7<r 



(53) 



Note A is a stochastic matrix and all diagonal entries sat- 
isfy An > 0. Hence the matrix is aperiodic which implies 
that the eigenvalue ax = 1 is the largest and is unique. 

We let Vi for i = 1, 2, . . . , N be the eigenvectors of A 
with corresponding eigenvalues |ai| > | a-2 1 > ■•■ > |ajv|- 
The eigenvector corresponding to the eigenvalue a\ = 1 
is Vi — (1,1,..., 1) T since the row-sums of A are equal 
to one. Recall that this means the distance vector dpc(n) 
does not tend to zero as n — > oo, but instead to a uniform 



phase shift ( 40 1 



lim dp C {s) =: 5 C 



(54) 



which has all components equal, (S ca ) i — 6^ for all i 
(i.e. all the neurons are at the same phase and hence in a 
globally synchronized state). Furthermore, recall that the 
distance from the globally synchronized state is given by 



dpc(nT) := max \S PC i(nT)\ 



(55) 



as defined in (411. Using the fact that ai = 1, V! = 

(1,1,..., 1) T and rewriting <5pc as a linear combination 
of the basis of eigenvectors gives 



N 



5'p C (nT) = 5p C (nT) - ^ = ^ fto?v< 



(56) 
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Then, since a 2 is the second largest eigenvalue, taking the 
infinity norm in ( 56 1 gives 



d PC (nT) 



max 

i 



\02 a 2 I max 

j 



V 2 



3 

N R 



\h\\a 2 \ 



max |v 2 j | 

j 



(57) 



where ~ means 'is asymptotically equal to (as n — > oo)' 
Taking the logarithm gives 



log (d PC (nT) ) - n log | a 2 \ 



(58) 



On the other hand dpc{n) asymptotically defines the syn- 
chronization time by 



dpc(nT) ~ exp - 



nT 
Tpc 



which after taking the logarithm gives 

nT 
Tpc 



log(d PC (n)) 



Comparing ( 58 1 with ( 60 1 leads to 

T 

TPC 



log |a 2 | 



(59) 



(60) 



(61) 



We numerically find the second largest eigenvalue of the 
matrix A and use this to calculate the analytical synchro- 
nization time ( |6f p . To minimize the influence of specific 
rewired networks and perturbations we average over 100 
realizations to obtain the average synchronization time 
(rpc). As for the Kuramoto system, the prediction of syn- 
chronization times based on the eigenvalues of the matrix 
A well agrees with those obtained from direct numerical 
simulation (Fig. [5} . 



5 The synchronization times for several 
network ensembles 

Here we study the dependence of the synchronization time 
on the topological randomness q and on the in-degree k. 
Following the original approach of Watts and Strogatz we 
first examine ensembles with a fixed in-degree k and then 
generalize to other generic network network ensembles. 
The results for different dynamics and different ensembles 
k(q) are depicted in Fig. [5] The synchronization times are 
studied with the emphasis on their behaviour in the small- 
world regime which is defined by Eq. ( 13 I and highlighted 
in Fig. [5] as the shaded region. 

Analytical predictions for the synchronization times in 
(19}, pi and ([611 are averaged over 100 network realiza- 



tions and depicted in Fig. [5] as solid lines. Simulation re- 
sults for the synchronization times are averaged over 100 
realizations of networks and perturbations. 



For the Kuramoto and pulse-coupled oscillators, de- 
termining the eigenvalues of the stability matrices of net- 
works yields synchronization time estimates that well agree 
with those found from direct numerical simulations. 

Only for the Rossler oscillators the synchronization 
times obtained from the numerical measurement of the de- 
caying maximal distances ^ show small but systematic 
deviations from the analytically predicted ones. These de- 
viations may be due to inaccuracies in determining the 
Euclidean distances that oscillate (see the decaying x- 
coordinates only in Fig.[l] b). 



5.1 Networks with fixed in-degree 

This subsection is dedicated to the dependence of average 
synchronization time (r (q, k)) on the topological random- 
ness q for standard Watts-Strogatz ensembles of networks 
with fixed in-degree k. 

We see in Fig. [l] and Fig. gb) that (t (q, k)) is mono- 
tonically decreasing with the topological randomness q 
and systematically depends on the network topology: Reg- 
ular ring networks (q — ¥ 0) are typically relatively slow to 
synchronize. We find that increasing q towards the small- 
world regime induces shorter and shorter network syn- 
chronization times, with small worlds synchronizing a few 
times faster than regular rings. Further increasing the ran- 
domness q induces even much faster synchronization, with 
fully random networks (q — > 1) synchronizing fastest (two 
orders of magnitude faster than small worlds in our ex- 
amples). Thus in network ensembles with fixed in-degree 
small worlds just occur intermediately during a monotonic 
increase of synchronization speed, but are not at all topo- 
logically optimal regarding their synchronization time. 

One could try to explain this dependence heuristically 
by the decrease of the average characteristic path length. 
Indeed, the dependence of (L (q, k)} on q is also monoton- 
ically decreasing in a similar fashion. It is intuitive that 
as the characteristic path length decreases, oscillators can 
communicate more efficiently and this leads to faster and 
more efficient synchronization. 



5.2 Ensembles with fixed average path length or 
betweenness centrality 

We therefore systematically study the synchronization time 
for generalized Watts-Strogatz ensembles of networks, spec- 
ified by a function k(q), where the average path length 
(L) is fixed while the degree of randomness q varies. We 
fix the average characteristic path length (L (q, k)) — 3.5 
as this gives us a wide range of q values. However, the 
results below are not sensitive to a change of this spe- 
cific value, cf. also [23 . We do not take k < 10 as the 
networks are in general no longer strongly connected for 
larger q values. For each of these in-degrees k a value of 
the topopolgical randomness q{k) is determined. Note that 
the standard deviations are larger for smaller q values. 
This is because we are rewiring a small number of edges 
here (Nk/q on average) and rewiring one edge more or 
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Fig. 5. (Color online) Comparing the synchronization times of different network ensembles, the small-world regime does not 
appear to be special, a: solid lines indicate network ensembles with fixed in-degree k = 50, with fixed average path length 
(L) — 3.5 and with fixed betweenness centrality (B) = 0.35. The small-world regime ( | 13| > is located between the red dashed 
lines, b - d: analytical predictions based on equations (18 1, (221 and (411 (solid lines) and simulation results (symbols with 



standard deviations) for synchronization times of Kuramoto (b 



ue, circle), Rossler (green, triangle) and pulse-coupled oscillators 



(red, square) for the sketched ensembles. 100 realizations were carried out in order to average over networks (and in simulations 
additionally over perturbations). 



less may have a strong effect on L(q 1 k) as it may add 
a long-range connection where there was not one previ- 
ously. Note that k decreases in a non-linear fashion as 
q increases for networks with (L (q, k)) = 3.5. When we 
increase q. we decrease the in-degree k. Thus, it might 
be expected that the amount of coupling each oscillator 
receives also decreases. This would affect the synchroniza- 
tion time [23 . We remove this factor by keeping the input 
each oscillator receives fixed ^ as q varies. By doing this, 
we have reduced the effect of changing the in-degree k on 
the synchronization time. Surprisingly, the synchroniza- 
tion time of network ensembles with fixed average char- 
acteristic path length non-monotonically depends on the 
topological randomness (Fig. [5] c). In particular, networks 
with intermediate randomness in the small-world regime 
synchronize slowest. 

Since faster synchronization times are apparently not 
only related to the decrease of the average path length, 
we investigated the dependence on other topological ob- 



servables which have been suggested to control whether 
or not a network actually synchronizes [13 , 43, 44, 45 , 46 
|4"7] , Representatively, an ensemble with fixed betweenness 
centralitjj^Jis shown in ((Fig. [5], d)). Ensembles with fixed 
clustering coefficient show a similar dependence, but also 
cover only a small range of q values. 



5.3 Synchronization time of generic network ensembles 

How does synchronization speed vary with randomness 
for more general ensembles k{q)l A systematic study of 
the synchronization time as a function of both in-degree k 
and randomness q (Fig. [6]) reveals an interesting nonlinear 
dependence. 

4 B(k, q) of node i (see e.g. |48j) measures the number of 
directed shortest paths that pass through this node, averaged 
over all nodes i. 
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Fig. 6. (Color online) Universal nonlinear dependence of syn- 
chronization time on in-degree k and topological randomness 
q. No generic ensemble k(q) exhibits fastest synchronization 
in the small-world regime, a: Kuramoto oscillators (modified 
from |22|). b: Rossler oscillators, c: pulse-coupled oscillators; 
logarithmic color scales. Parameters chosen as in Fig.JTj Syn- 
chronization times are obtained from equations (19 1, (31 1 and 

dsn. 



Firstly, it confirms that for all networks with fixed in- 
degree k the synchronization time is monotonic in the 
randomness q and the small-world regime at intermedi- 
ate randomness is not specifically distinguished. 

Secondly, the two-dimensional function (r(q, k)) im- 
plies that ensembles of networks with fixed path lengths 
all exhibit a non-monotonic behavior of the synchroniza- 
tion time, with slowest synchronization for intermediate 
randomness. 

Thirdly, considering graph ensembles characterized by 
any other smooth function k(q), q £ [0,1], shows that 
this is a general phenomenon and the specific choice of an 
ensemble k(q) is not essential. 



In fact, as illustrated in Fig.[6j for any generic network 
ensemble k(q) (including ensembles with fixed in-degree, 
fixed path length and fixed betweenness centrality as spe- 
cial choice) the synchronization speed (r(q, k(q))) is either 
intermediate or slowest, but never fastest at intermediate 
randomness, in particular in the small- world regime. It 
is remarkable that this seems to hold universally as the 
synchronization times are similar for Kuramoto oscillators 
(Fig. [6] a), periodic (Fig. [6 b) and chaotic Rossler 
oscillators and pulse-coupled oscillators (Fig. pi c) . 



5.4 Relation between Kuramoto and pulse-coupled 
oscillators 

Comparing the synchronization times for Kuramoto and 
pulse-coupled oscillators in Fig. [6] both show a striking 
similarity: Could the synchronization times be mapped 
on each other? 

Therefore let us investigate how perturbations 6(t) 
evolve in both systems: For the Kuramoto oscillators we 
have 



S K (t) 



„Rc \ 2 t 



5(0) 



(62) 



while perturbations in the pulse-coupled system propagate 
like 

S PC (nT) = \a 2 \ n S(0). (63) 

Obtaining similar synchronization times for both dynam- 
ics demands these perturbations to evolve in the same way. 
Setting nT := t the crucial eigenvalues should satisfy 



Re A 2 = 



log i j 
T 



(64) 



Comparing the structure of the two relevant matrices A 
( 52 1 and A pi , we obtain the following relation for the 



respective eigenvalues a 2 and A2 

ck 



A, 



cpc 



a 2 - 1 



(65) 



where the quotient ck / cpc depends on the system param- 
eters 



ck =o"k 
cpc = 



7°pc 
Ie-y A + 7CTpc ' 



(66) 
(67) 



Note that <tk > while trpc < 0. Setting the parameters 
in the way that this quotient and the period T (T is close 
to one anyhow) are both equal to one, we get 

1 + log |a.2 1 = Rea2 provided \a 2 \ « 1 , Rea2 ~ 1. (68) 

This means that synchronization times obtained for net- 
works of pulse-coupled oscillators and for the same struc- 
tured networks of Kuramoto oscillators are equal if the 
second largest eigenvalue of matrix A satisfies Re a 2 ~ 1 
and \a 2 \ ~ 1, in the sense that |log|a 2 || <C 1. This is in 
general the case for large networks, for which the stochas- 
tic matrix A has TV eigenvalues with real parts in [—1,1]. 
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(air transport network: the 
network: C. Elegans |51| . or§ 
on the trajectories of dollar 
d: synchronization times (t) 
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-world networks consistently synchronize several times slower than their randomized counterparts 
US airports with the largest amount of traffic |49j . social network: inmates in prison [50 , neural 
;anisational network: research team in a manufacturing company |52| . human travel network: based 
bills [53]). a: clustering coefficients (C). b: average path lengths (L). c: synchronization times (t) 
relative to the average path length (L). Inset: the synchronization times for C. Elegans divided by 
show slightly non-monotonic behaviour. All quantities are averaged over 100 network realizations. 



6 Real- world networks 

For Watts-Strogatz small-world model networks we have 
found that the synchronization speed is either intermedi- 
ate or slowest, but never fastest in the small- world regime. 
Moreover, keeping the in-degree fixed, the model networks 
synchronize the faster the more random they are. To sup- 
port that this monotone relation also holds in much more 
generality we considered various real-world networks: an 
air transport network a social network [5U], a neural 
network |51j . an organisational network [52J and a human 
travel network |53| . 

Randomizing is performed by rewiring the existing con- 
nections with rewiring probability p G [0, 1]. Note that the 
rewiring probability p and the topological randomness q 
are two different quantities. Here we start with the original 
real- world network (p = 0), which may be in the small- 
world regime already. The rewiring process is performed 
as explained in Section [3] only the outgoing edges of the 
directed links are rewired, which means that each node's 
in-degree stays constant during the rewiring. Having con- 
sidered networks with the same in-degrees for all nodes 
so far, these in-degrees may vary in real-world networks: 
But due to the rewiring routine the network's initial in- 



degree distribution is kept constant during rewiring. Thus 
it is not uncommon that a network may split from only 
one strongly-connected component (SCC) in the begin- 
ning to several ones |25j . If splitting occurs the rewiring is 
repeated until the resulting network consists of one SCC 
only again. All measured quantities are averaged over 100 
network realizations. 

The studied real-world networks show small-world be- 
haviour: gradually randomizing these networks, their clus- 
tering coefficients and their average path lengths mono- 
tonically decrease (Fig. [7] a and b). Only the cluster- 
ing coefficient of the human travel network partly shows 
a non-monotonic behaviour. This network is extremely 
large (464670 nodes) in comparison to the other networks 
(w 100 - 500 nodes). 

We found that all real-world networks consistently syn- 
chronized several times slower than their randomized coun- 
terparts (Fig. [7] c). For all networks, the synchroniza- 
tion times monotonically increase with increasing random 
rewiring. For all but the neural circuit of C. Elegans |51| . 
this holds even for ensembles with fixed average path length; 
for the latter we observed non-monotonic dependence with 
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slowest synchronization for intermediate randomness (Fig.|7| 
d). 



7 Conclusion 

We investigated the effect of topology on the speed of syn- 
chronization of various oscillator types, intrinsic dynamics 
and coupling schemes: phase oscillators coupled via phase 
differences, higher-dimensional periodic systems coupled 
diffusively as well as neural circuits with inhibitory de- 
layed pulse-coupling and consistently found qualitatively 
the same results. We derived analytical predictions for the 
asymptotic synchronization times, including an extension 
of the master stability function to determine how fast the 
system actually synchronizes. We compared the synchro- 
nization speed for different network ensembles: 

Firstly, we found that for networks of fixed in-degree 
fc, the average synchronization time (r) is monotonically 
decreasing with the topological randomness q. Comparing 
different fixed- fc ensembles small- world networks always 
synchronize quicker than regular networks. 

Secondly, the intuitive idea that this is due to the 
decrease in average characteristic path length (L(q, fc)) 
does not provide a complete explanation: Instead of fix- 
ing the in-degree, we fixed the average characteristic path 
length. For such ensembles networks in the small-world 
regime synchronize slower than regular networks and the 
synchronization speed non-monotonically depends on the 
topological randomness q. The in-degree fc is monotoni- 
cally decreasing with q and so does the clustering coeffi- 
cient (C(q, fc)). So neither of these topological properties 
alone gives an obvious explanation and the phenomenon 
results from an interplay between several network prop- 
erties. For example, the faster synchronization of regular 
networks than of small-world networks may be due to the 
in-degree fc being large. This is not because the total cou- 
pling strength = fccr is high, as we kept this fixed 
for all (fc, g)-pairs, but may be because the oscillators re- 
ceive the coupling effect from a large number of oscillators. 
However, we also see fast synchronization for random net- 
works where fc is small and so the same total coupling 
amount is received from far fewer interacting oscillators. 
So the explanation for the non-monotonic dependence is 
non-trivial. 

Thirdly, we investigated the dependence on other topo- 
logical observables apart from small-world properties: net- 
work ensembles with fixed betweenness centrality have 
been displayed as an example, but yet a simple expla- 
nation for the nonlinear dependence is missing. 

Further, we studied the full nonlinear dependence of 
the synchronization time on the in-degree fc and the topo- 
logical randomness q for generic network ensembles. We 
found that fastest synchronization is essentially impossible 
in the small- world regime, except for highly artificial en- 
sembles. This statement holds for all observed dynamics. 
In particular, the synchronization times for the Kuramoto 
and pulse-coupled oscillators are strikingly similar. 

ft would be interesting to extend the analysis started 



chronization times for Kuramoto oscillators could be ap- 
proximated by or even analytically mapped on the times 
for the pulse-coupled oscillators. Additionally an under- 
standing or even an analytical description of the curves 
of same synchronization times in Fig. [6] is extremely help- 
ful for finding further relations between the topology and 
the dynamics of complex networks. In particular, it is an 
interesting question to understand the behaviour of the 
second largest eigenvalue of the Laplacian as a function of 
q for fixed fc and N. First results indicate that this might 
follow a simple power law with exponent close to 1 for 
values of q > 1/N. 

In this article we focused on systems with fixed size N. 
In the Watts-Strogatz ensemble the scaling of quantities 
such as the average path length L or the clustering coeffi- 
cient C with the system size depends heavily on q [6], e.g. 
L ~ N for q — and L ~ log TV for q — 1. Therefore, it 
would be an interesting question to study how the results 
illustrated in Fig. [6] change with N and what would be 
the appropriate definition of the small-world regime and 
other generalized ensembles with given structural features 
as a function of the system size. 

We found that small worlds in general never synchro- 
nize fastest. Specifically, in networks with fixed average 
path length, synchrony is consistently fast for regular rings, 
fastest for completely random networks, and slowest in the 
intermediate small- world regime (Fig. [SJ. It is an aston- 
ishing result that this holds across various oscillator types, 
intrinsic dynamics and coupling schemes: phase oscillators 
coupled via phase differences, higher-dimensional periodic 
systems coupled diffusively as well as neural circuits with 
inhibitory delayed pulse-coupling. 

In particular, small- world topologies are not at all spe- 
cial and may synchronize orders of magnitude slower than 
completely random networks. So generically the small- 
world regime either exhibits slowest synchronization or 
just exhibits no extremal properties regarding synchro- 
nization times. 

Given the variety of the investigated dynamical sys- 
tems our results indicate that this is a universal phe- 
nomenon. 

Our investigations of real-world networks support this 
view. Although the considered networks may be in the 
small-world regime already, rewiring still strongly increases 
the synchronization speed, even for ensembles with fixed 
average path length. It remains an open question why 
rewiring typically implies faster synchronization. 

We thank C. Kirst and D. Brockmann for providing 
data of real-world networks. Supported by the BMBF un- 
der grant no. 01GQ1005B, by a grant of the Max Planck 
Society to MT and by EPSRC grant no. EP/E501311/1 to 
SG. SG is also grateful for the hospitality of the Hausdorff 
Research Institute for Mathematics in Bonn. 



in Section 5.4 to find out under which conditions the syn- 
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